Diffusion

import numpy as np
import scipy 
import ipywidgets as widgets
import matplotlib.pyplot as plt

%matplotlib inline

Microscopic aspects of Diffusion

• Einstein developed a theory of diffusion based on random walk ideas and obtained a key equation relating mean square displacement of particle in d dimensions to time \((t = n \Delta t)\). Generate N random walks with \(n\) steps all starting from origin \(r_0=0\).

\[R_n = \sum^{n}_{i=0} r_i\]

• Expressing steps in terms of time increments \(n=\frac{t}{\delta t}\) we compute average over N random walkers in \(d=1,2,3\) dimenion.

\[\langle R^2_n\rangle = \sum_i \sum_j \langle r_i r_j \rangle = \sum_i \langle r^2_i \rangle = \sum_i d \cdot \langle \delta x^2_i \rangle = d \cdot n \cdot \delta x^2\]

\[\langle R^2_n\rangle = d \cdot \frac{t}{\delta t} \cdot \delta x^2\]

• Grouping constants together we define the diffusion coefficient which is expressed in terms of microscopic quantities defined in the random walk model!

\[\langle R^2 (t) \rangle = 2d \cdot D \cdot t\]

\[D = \frac{\langle \delta x^2 \rangle}{2\delta t}\]

• We end up with a general expression for a mean square displacement as a function of time. Any motion which adheres to this scaling with time will be called diffusive.

\[MSD(t) = 2d D \cdot t^{1/2}\]

d=3 
t= np.linspace(0, 1, 10000)

for D in [0.01, 0.1, 1, 10]:
    
    plt.plot(t, 2*d *D * t**0.5)
    #plt.loglog(t, 2*d *D * t**0.5)

plt.ylabel('$MSD(t)$')
plt.xlabel('$t$')

Text(0.5, 0, '$t$')

Brownian motion

• The Brownian motion describes the motion of a particle suspended in a solvent consisting of much smaller moleules. The displacement of particle is generatd by a sum of large number of independent random collisions with solvent molecules. Hence we invoke Central Limit Theorem to approximate displacement at each \(dt\) as a normalally idstributed random variable:

\[x(t+dt)-x(t)=N(0,\sqrt{2D dt})\]

• We assume we have started at position \(\mu=0\), and our variance is given by \(\sigma^2=2Dt\), Where D is the diffusion coefficient, which is related to the parameters of the discrete random walk as shown in the lecture.

\[x(t+dt)=x(t)+\sqrt{2D dt} \cdot N(0,1)\]

• In the last step, we re-wrote Brownian motion equation in a convenient way by shifting the normally distributed random variable by mean and scaling by standard deviation \(N(\mu, \sigma^2) = \mu + \sigma N(0,1)\)

def brown(T, N, dt=1, D=1, d=3):
    
    """
    Creates 3D brownian path given:
    time T 
    N=1 trajecotires
    dt=1 timestep
    D=1 diffusion coeff
    returns np.array with shape (N, T, 3)
    """
    
    n = int(T/dt) # how many points to sample
    
    dR = np.sqrt(2*D*dt) * np.random.randn(N, n, d) # 3D position of brownian particle
    
    R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
    
    return R

R=brown(T=3000, N=1000)
print(R.shape)

(1000, 3000, 3)

def brownian_plot(t=10):
    
    fig, ax = plt.subplots(ncols=2)
    
    ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
    
    ax[1].hist(R[:, 10, 0], density=True, color='red')
    ax[1].hist(R[:, t, 0], density=True)
    
    ax[1].set_ylim([0,0.1])
    
    ax[0].set_ylim([-200, 200])
    ax[0].set_xlim([-200, 200])
    
    fig.tight_layout()

import holoviews as hv
hv.extension('plotly')

plots = []
for i in range(10):
    
    plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
    plots.append(plot)
    
hv.Overlay(plots) 

rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]

xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])

hv.Overlay(rw_curve) << ydist << xdist

Diffusion Equation

The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\). Formulated empirically as Fick’s laws

\[\frac{\partial\rho}{\partial t} = \mathcal{D}\nabla^2\rho\]

• This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour

• This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.

• Can solve numerically by writing derivatives as finite differences!

• Can also simulate via random walk!

• Diffusion coefficient \(D\), Units \([L^2]/[T]\)

Important special case solution (here written in 1d):

\[\rho(x,t) = \frac{1}{\sqrt{2\pi \sigma(t)^2}}\exp\left(-\frac{x^2}{2\sigma(t)^2}\right),\]

where \(\sigma(t)=\sqrt{2{D}t}\)

• density remains Gaussian and Gaussian becomes wider with time1

• check that this is indeed a solution by plugging into the diffusion equation!

def sigma(t, D = 1):
    return np.sqrt(2*D*t)

def gaussian(x, t):
    return  1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #

def diffusion(t=0.001):
    
    R = brown(T=101, N=1000)
    x = np.linspace(-20, 20, 100)
    
    plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
    
    plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
    
    plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
    
    plt.legend()
    plt.ylabel('$p(x)$')
    plt.xlabel('$x$')
    plt.xlim([-25, 25])